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Abstract. We numerically investigate Lyapunov instabilities for one-, two- and 
three-dimensional lattices of interacting classical spins at infinite temperature. 
We obtain the largest Lyapunov exponents for a very large variety of nearest- 
neighbor spin-spin interactions and complete Lyapunov spectra in a few selected 
cases. All Lyapunov spectra have weakly convex shapes and do not exhibit the 
Lyapunov modes observed for gases of hard-core particles. We investigate the 
dependence of the largest Lyapunov exponents and whole Lyapunov spectra on 
the lattice size and find that both quickly become size-independent. Finally, we 
analyze the dependence of the largest Lyapunov exponents on the anisotropy of 
spin-spin interaction with the particular focus on the difference between bipartite 
and nonbipartite lattices. 
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1. Introduction 

Investigations of the Lyapunov instabilities in many-particle systems are often 
motivated by the role of chaos in the foundations of statistical physics [1 , % 3 , which is 
still not fully understood. In general, interacting, many-particle, classical systems are 
expected to be chaotic. The largest Lyapunov exponents and properties of the entire 
Lyapunov spectra have been calculated numerically for classical many-body systems, 
such as gases of hard-core particles [U [5] , fluids with soft interactions [5J , and lattice 
two-dimensional rotators [BJ [7] , and analytically in a few cases [HI M Q33 H3 HH EH H3 ■ 
In the present paper, we focus on lattices of interacting classical spins. 

Classical spins often appear in the theoretical studies as the large-spin limit of 
quantum spins. Moreover, even when one deals with lattices of spins 1/2, parallels 
between classical and quantum spin dynamics still remain. These parallels have 
recently received much attention, in particular, in the context of the asymptotic 
exponential-oscillatory behavior of nuclear spin decays in solids [HI [HI [TT1 [HI [TH 
[201 121] • These decays was identified in Refs. [HUH] with chaotic eigenmodes in both 
classical and quantum many-spin systems. 

Although it appears very likely a priori that lattices of interacting classical spins 
exhibit chaotic dynamics, no systematic investigation of the chaotic properties of these 
lattice was undertaken until our previous work [22) , which presented a survey of the 
largest Lyapunov exponents for a very large variety of spin lattices and Hamiltonian 
anisotropies. The principal finding of Ref. [22] was that all Hamiltonians considered, 
with the exception of the Ising case, led to chaotic dynamics as evidenced by the 
nonzero value of the largest Lyapunov exponent. We also obtained both analytically 
and numerically the power- law scaling of the largest Lyapunov exponent in the vicinity 
of the integrable Ising limit. 

In the present paper, we complement the findings of Ref. [22] in several 
respects. Namely, we compute complete Lyapunov spectra for a few selected 
spin lattices and show their dependence on the lattice size, and also present a 
more extensive investigation of the lattice size dependence of the largest Lyapunov 
exponent. We discuss the dependence of the largest Lyapunov exponent on the 
Hamiltonian anisotropy with particular emphasis on the difference between bipartite 
and nonbipartite lattices. 

2. General formulation 

2.1. Spin model 

We consider periodically closed spin lattices with the nearest-neighbor (NN) 
interaction Hamiltonian of the following kind: 

XN 



where (Si X , Si y , Si Z ) = Si are the three projections of the classical spin vector 
of unit length on the ith lattice (i.e. = 1), and J x , J y , J z are the coupling 
constants, which we also normalize by condition J% + Jy + J? = 1- Below, we 
often mention Ising, Heisenberg and "anti-Heisenberg" limits of the Hamiltonian 
(JTJ. The Ising Hamiltonian corresponds to 4 = J y = 0, J z = 1, Heisenberg 
Hamiltonian to J x = J y = J z = 1/V3, and, finally, anti-Heisenberg Hamiltonian 
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Figure 1. Lattices investigated in this work. Bipartite: (Ll) chain, (L2) 
rectangular ladder, (L3) square lattice, (L4) bilayer of square lattices, (L5) cubic 
lattice. Non-bipartite: (L6) triangular ladder, (L7) triangular lattice. 



to J x = J y = — J z = 1/V3- The total number of spins in the lattice is denoted as N. 
The phase space of such a lattice has dimensionality 27V. 

We consider seven lattices shown in Fig. Q] and labeled as (L1-L7). Lattices (Ll- 
L5) are bipartite, which means that they can be divided into two sublattices such that 
all the interacting neighbors for a spin on one sublattice belong to the other sublattice. 
Lattices (L6,L7) are nonbipartite. 

The equations of motion associated with the Hamiltonian ([1]) can be obtained in 
the Poisson-bracket formalism: dSi^/dt = {%, Si,^}, where the index \i admits values 
1, 2 or 3 representing the projections x, y, or z, respectively. The primary Poisson 
brackets are: {Si ifl , Sj M } = Sij J2 K e ni/K^i,Ki where Sij is the Kronecker delta and e^^ 
the Levi-Civita symbol. The resulting equations of motion thus become 

Si = Si x h, , (2) 

where is the local field given by expression 

h-i ^ ^ JxSjx^x H~ Jy^ 1> jy^y J zS j z^z ■ (3) 

M 

Here e x , e y and e z are the unit vectors along the respective directions, and j(i) implies 
the summation of the nearest neighbors of the i-th. lattice site. 

In this work, we restrict ourselves to the Lyapunov instabilities on the zero energy 
shell, which corresponds to infinite temperature in the microcanonical sense. 



2.2. Analytical considerations 

Since the dynamics of this system are fully time-reversible, positive and negative 
Lyapunov exponents are expected to form conjugate pairs with equal absolute values. 
In the general case of different J x , J y , and J z , there should be only two zero Lyapunov 
exponents corresponding to the energy and the time shift directions. In the case 
Jx — Jy 7^ Jz, of which the anti-Heisenberg case is an example, the z-component of 
the total spin polarization becomes the integral of motion, and hence one more pair 
of Lyapunov exponents should also become equal to zero. In the Heisenberg case, 
Jx = Jy = Jzi all three components of the total spin polarization become integrals of 
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motion. However, these three integrals of motion are not dynamically independent, 
because one of them can be obtained as the Poisson bracket of the two others. The 
two independent integrals of motion thus imply two additional pairs of zero Lyapunov 
exponents. In the Ising case, the z-component of each spin is the integral of motion. 
Hence, the system is integrable and all Lyapunov exponents are equal to zero. 

Now we discuss the connections between the Lyapunov spectra of different 
anisotropic Hamiltonians. We first note that infinite-temperature Lyapunov spectra 
are expected to be identical for the Hamiltonians with the coupling constants 
(J x ,J y ,J z ) and (—J x ,—J y ,—J z ), because the zero-energy shells in the two cases 
are identical, while the change of the sign of the coupling constants amounts to 
the operation of time-reversal. We further note that, for bipartite lattices, the 
infinite temperature Lyapunov spectra for the coupling constants (J x , J y , J z ) and 
(—Jx,—Jy,Jz) should also be identical. In order to see this, one should make the 
transformation Six ~^ ~Si X and Si y — > — Si y for one of the two sublattices forming 
the bipartite lattice and then examine the resulting equations of motion. 

The two observations made above imply that, for a given bipartite lattice, 
the Lyapunov spectra for the Heisenberg and the anti-Heisenberg Hamiltonians are 
identical to each other, which means that, in the latter case, the Lyapunov spectrum 
has three pairs of zero Lyapunov exponents instead of two expected for the generic case 
of J x = J y =/= J z . The extra pair of the zero exponents originates from the following 
two (not independent) integrals of motion: rfSix or J2i(~l)^Si y , where £ is the 

index taking value for one sublattice and 1 for the other. The same considerations 
also imply that any fully anisotropic Hamiltonian with J x = — J z or equivalent on a 
bipartite lattice can be converted to the axially symmetric Hamiltonian with J x = J z . 
Therefore, the original Hamiltonain has an extra pair of zero Lyapunov exponents 
corresponding to the integral of motion ^2,A— 1) ^Si y . 

Small spin clusters may have additional nontrivial integrals of motion. In partic- 
ular, any 4-spin periodic chain with the general anisotropic Hamiltonian of form (flj is 
fully integrable. This is because the first and third spins in this chain rotate in the same 
local field [J x (S2x + S^), J y (S2 y + S^ y ), J Z (S2 Z + S^)], while the second and the 
fourth spins rotate in another local field [J x (Si x + S^), J y (S\ y + Ss y ), J Z (S\ Z + S^ z )]. 
Therefore, there arc two additional integrals of motion, namely: Si • S3 and S2 • S4. 
Finally, one can also check that (Si + S3) • (S2 + S4) is also the integral of motion. As a 
result, the number of the integrals of motion (including energy) becomes 4, while the 
dimensionality of the phase space is 8, i.e. the problem is fully integrable. In the case 
of the isotropic Heisenberg Hamiltonian, a much larger variety of small spin clusters 
with nontrivial integrals of motion were cataloged in Ref. [23] . 

2.3. Numerical simulations 

The equations of motion were integrated using the fourth-order Runge-Kutta 
algorithm with a time-step St = 0.005. This is sufficiently small so that on the time 
scale of our simulations, energy is conserved with the 6-digit accuracy. Simulations 
were run for 20000 time units, which was sufficient for accurate convergence of the 
smaller exponents. 

The Lyapunov exponents were obtained using a version of the standard 
reorthonormalization algorithm [4, 24, 25]. In order to obtain first n largest Lyapunov 
exponents, we numerically propagate a reference trajectory f(t) and n initially 
orthogonal perturbation vectors S-y^t). At every time step, we propagate the 
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perturbations using the linear tangent space map that is obtained by numerically 
taking derivatives, 

8j(t + to)= *Mp^ 5>rW. (4) 

After each time interval At — 0.25, we hierarchically reorthogonolize the perturbation 
vectors 5'y i (t) using the Gram-Schmidt procedure, and then renormalize their lengths 
back to the initial values. The renormahzation factors are recorded for the kth time 
interval At as ati(k). Finally, we compute the ith Lyapunov exponent using the formula 
A,; = -j^i Sfc=i lnat(^), where K is the number of the renormahzation time intervals. 

As the test of the accuracy of our numerical routine, we have checked that the 
Lyapunov exponents form conjugate pairs and that the number of the zero Lyapunov 
exponents is equal to the number expected from the symmetry of the Hamiltonian as 
discussed in Section 12.21 We have also checked the symmetries between Heisenberg 
and anti-Heisenberg Lyapunov spectra expected based on the arguments presented in 
Section O 

In order to obtain initial conditions at zero total energy, we first choose random 
orientations for all spins. This produces a total energy close to zero, but with 
fluctuations of the order of y/~N. Then, in order to arrive to the zero total energy, we 
evolve the dissipative equations of motion: 

Si = ±Si x (Si x hi) . (5) 

Depending on the sign in front of the right-hand side, this increases or decreases the 
total energy associated with the Hamiltonian ([T]). Once the zero value of the total 
energy is reached, we additionally assure the randomness of the initial conditions on 
the energy shell by performing I0A sequential rotations of random spins by random 
angles around the directions of their respective local fields hi given by Eq.Q. These 
rotations preserve the total energy. 



3. Results and discussion 



3.1. Lyapunov spectra 

We have computed the full Lyapunov spectra for the linear chain (LI), cubic lattice 
(L5) and triangular lattice (L7) with the Heisenberg Hamiltoninan, and also for 
the triangular lattice (L7) with the anti-Heisenberg Hamiltonian. The results are 
presented in Fig. [2] The spectra are typically weakly convex, regardless of the lattice. 
In the case of the cubic lattice with Heisenberg interaction, the spectrum is actually 
very close to linear. As explained in Section l2"T2l the Lyapunov spectra for the bipartite 
lattices (LI) and (L5) with anti-Heisenberg Hamiltonian are identical to those already 
presented in Fig. [2] for the Heisenberg Hamiltonian. The effect of the lattice size on 
the Lyapunov spectrum for the linear chain with the Heisenberg Hamiltonian is shown 
in Fig. |3l The size dependence of the spectrum appears to be very weak for N > 4 
and becomes virtually unobservable for the chains containing more than 32 spins. 

Comparing the Lyapunov spectra of classical spins with the Lyapunov spectra 
of other many-particle systems, we first note, that the spectra presented in Fig. [2] do 
not exhibit an offset from zero for the smallest positive exponents. Such an offset 
was observed in gases of particles [4] and high-dimensional billiards [26l H2] but not 
in systems with sufficiently soft interactions [B]. The classical spin lattices obviously 
belong to the latter group. We further note that delocalised Lyapunov-Goldstone 
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Figure 2. (Color online) Lyapunov spectra for lattices (LI), (L5), and (L7) 
with Heisenberg and anti-Heisenberg Hamiltonians. The number of spins in each 
case is N = 256. Lyapunov exponents A; are ordered according their values and 
plotted as a function of i: (a) Positive halves of the spectra. The inset shows 
the same plots but with the vertical axis renormalized by the value of the largest 
Lyapunov exponent Ai . (b) Magnified small-exponent parts of the spectra. Note: 
lattices (LI) and (L5) have identical Lyapunov spectra for the Heisenberg and 
the anti-Heisenberg interactions (see Section 12.21 1. In (b), the numerical error is 
roughly the same as the size of the symbols. 




Figure 3. (Color online) The positive halves of the Lyapunov spectra for spin 
chains (LI) of different lengths TV with the Heisenberg Hamiltonian. Note: for 
N = 4, all Lyapunov exponents are zero, due to the symmetries in the system 
(see Section [2~2t . 



modes, that were observed in dilute gases [U [57J [25J [TT] and in some other extended 
systems [29], were not observed in the present case. 

3.2. Largest Lyapunov exponents: finite size effects 

Accurate numerical calculations of full Lyapunov spectra are very demanding, because 
the computational cost grows as N 2 . For this reason, the lattices investigated in 
the preceding subsection were relatively small. In this subsection, we focus only on 
the largest Lyapunov exponents, Ai. The cost of computing Ai grows only as N, 
which allows us to investigate much larger lattices and a much greater variety of the 
Hamiltonians. An extensive investigation of this kind was already reported by us in 
Ref. [22] . In the present and the next subsections we present some results and analysis 
that were not included in Ref. [22] . 
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Figure 4. (Color online) The dependence of the largest Lyapunov exponent on 
the number of spins N for each of the seven lattices indicated in the plot legent for 
(a) Heisenberg coupling and (b) anti-Heisenberg coupling. For bipartite lattices, 
the largest Lyapunov exponents are the same for Heisenberg and anti-Heisenberg 
coupling (see Section l2.2t . The largest Lyapunov exponents become constant for 
sufficiently large systems. 



In this subsection, we investigate the dependence of Ai on the lattice size for all 
seven lattices shown in Fig. [1] with Heisenberg and the anti-Heisenberg Hamiltonians. 
The results are shown in Fig. 01 They indicate that, within our numerical accuracy, Ai 
becomes size-independent for sufficiently large lattices. In particular, on the basis of 
these results even slow logarithmic growth of Ai with the lattice size can be excluded. 

The saturation of Ai with the lattice size can be explained on the basis of the 
following consideration. The exponential growth of the perturbation vector 5~y 1 in 
many-spin phase space with rate Ai implies that the projection Sf 1 on the subspace 
of each individual spin {Si x , Si y , Si Z } should also, on average, grow exponentially with 
the same rate. However, the instantaneous growth rate for the perturbation of the 
coordinates of a given spin is limited from above by the value of the order of the 
maximum value of the local local field max(|h,-|) = no max(| J x \, \J y \, \J Z \), where no 
is the number of nearest neighbors for each lattice site. Since this constraint does not 
depend on the lattice size, the growth of the largest Lyapunov exponent must saturate 
as the lattice size increases. 

3.3. Largest Lyapunov exponents: dependence on the Hamiltonian anisotropy 

In Ref. [22], we conducted a systematic survey of the dependence of Ai on the 
Hamiltonian anisotropy on the "interaction sphere" constrained by the condition 
J\ ' + Jy + = 1. We have found that the principal parameter controlling this 
dependence is J max = max(| J x \, \J y \,\J z \). In particular, this parameter quantifies the 
approach to the integrable Ising case corresponding to J max = 1. The main plot of 
Ref. [22] is reproduced in Fig. [5] 

Here we focus on the difference between the anisotropy dependence of Ai for 
bipartite lattices (L1-L5) and nonbipartite lattices (L6, L7). As can be seen in Fig.[S] 
the bipartite lattices (L1-L5) show a nearly universal dependence Ai(J max ), which 
scales only with the number of interacting neighbors. The small spread of the sampled 
values of Ai at a given value of J max indicates that Ai depends only very little on the 
ratio of the two coupling constants that have smaller absolute values. Our investigation 
of the nonbipartite lattices (L6,L7) was, in fact, motivated by the impression that the 
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Figure 5. (Color online) Largest Lyapunov exponents as presented in Ref. |22j . 
Each point represents one Ai obtained numerically for a lattice indicated in the 
plot legend with one randomly chosen set of values J x , J y and J z as described in 
the text. 

number of interacting neighbors alone determines the entire anisotropy dependence of 
Ai. We wanted to compare Ai for the lattices (L6) and (L3), where, in the both cases, 
each site has four nearest neighbors, and for the lattices (L7) and (L5), each site has 
six nearest neighbors. 

We found, however, that, as seen in Fig.[5j the nonbipartite lattices (L6) and (L7) 
exhibit a noticeable fork-like spread of Ai as Jmax approaches l/y/3. The upper and the 
lower tips of the fork correspond to the anti-Heisenberg and Heisenberg Hamiltonians, 
respectively. In the anti-Heisenberg case, the value of Ai is close to that of a bipartite 
lattice with the same number of nearest neighbors. In the Heisenberg case, the value 
of Xi is closer to the bipartite lattice with one fewer nearest neighbor. This spread 
indicates that the knowledge of J max alone is insufficient to determine Ai. The two- 
dimensional anisotropy dependence behind this spread is shown in Fig. [SI where we 
present it for lattices (L5) and (L7) in the form of the color density plots as a function 
of J x and J y . 

Less obvious from Fig. [5j is the fact that the significant majority of the sampled 
values of Ai for the lattices (L6) and (L7) agree very well with the dependence 
Ai(Jmax) for the bipartite lattices (L3) and (L5), respectively. The comparison 
between Figs.|B(a) and (b) clearly shows that the difference between lattices (L5) and 
(L7) is only pronounced, when all three interaction constants have the same sign and 
roughly the same value — or, in other words, when they approach the Heisenberg limit. 
This implies that Ai for the anti-Heisenberg Hamiltonian has a more typical value than 
for the Heisenberg Hamiltonian, and that the expected universality of Ai( J max ) for the 
same number of nearest neighbors still roughly holds even for nonbipartite lattices. 
It thus appears that the conservation of the total value of Sj in the Heisenberg case 
leads to the reduction of the effective number of the nearest neighbors by roughly one, 
as far as the value of Ai is concerned. 

We suspect that the situation here is similar to the origin of the frustrated 
low-temperature magnetism for the Heisenberg model on nonbipartite lattices. The 
interaction energy for a spin pair is minimal when the two spins are antiparallel to 
each other, but, on a nonbipartite lattice, such an antiparallel configuration cannot 
simultaneously exist for all pairs of interacting spins. Hence the ground state of 
a nonbipartite lattice is frustrated. In the case of the Lyapunov instabilities, the 
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Figure 6. (Color online) Largest Lyapunov exponent as a function of the 
coupling constants J x and J y . The third coupling constant is J z = y 1 — Jj — J ^ . 
Plots (a) and (b) represent, respectively, the numerically generated results for 
lattices (L5) and (L7). Both lattices have six nearest neighbors, but lattice (L5) 
is bipartite, while lattice (L7) is not. As explained in Section 12.21 the plots for 

J z = — yjl — J% — J2 can be obtained from the plots shown in this figure by 
reversing the signs of J x and J y . 

conservation of the total spin polarization implies that the perturbation vector 5'y 1 
corresponding to the largest Lyapunov exponent does not grow along the direction 
of the total spin polarization of the entire lattice. Locally this means that, if the 
projection of 6~f 1 on a given spin's subspace grows, then it should be compensated by 
the growth in the opposite directions for the projections of Sj-^ onto the subspaces of 
the nearest neighbors. Achieving the maximum growth of the perturbation, and hence 
the largest value of the Lyapunov exponent presumably requires that the perturbation 
vector 6~/i maximizes this antialignment of the projections onto adjacent sites. For 
the bipartite lattices, the antialignment can, in principle, be made perfect, but for 
the nonbipartite lattices this is impossible. Such an explanation is consistent with the 
fact that Ai for lattices (L6,L7) in the Heisenberg limit is smaller than in the anti- 
Heisenberg limit. It seems also to be connected to the fact that the above reduction 
approximately leads to the value of Ai for a bipartite lattice with roughly one fewer 
nearest neighbor per site. 

We finally remark that the overall small spread of values of Ai for bipartite 
lattices at a given value of J ma x is related to the symmetries described in 
Section 12 .21 which imply that the Lyapunov exponents are identical for eight 
combinations of the coupling constants characterized by the same value of J ma x, 

namely. (Jx: Jy: Jzjt ( Jx: Jy:Jz): ( JxtJy: Jz): (Jx: Jy: Jz): ( Jx: Jy: Jz): 

(—Jx:Jy:Jz): {Jx:~Jy: Jz): (Jx, Jy: —Jz)- The values of Ai cannot change much 
between these eight points. 

4. Summary and conclusions 

We have investigated numerically the Lyapunov spectra of systems of many classical 
spins for a variety of lattices and coupling constants at infinite temperature. The 
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possibility of varying the coupling constants, from the highly symmetric isotropic 
Heisenberg model, through partially symmetric couplings such as the anti-Heisenberg 
model, to the completely anisotropic integrable case, makes these systems particularly 
interesting. We have presented: (i) calculations of the Lyapunov spectra for selected 
lattices of interacting classical spins; (ii) investigations of the lattice-size dependence; 
(iii) investigations the largest Lyapunov exponents for a broader group of coupling 
constants, lattices and large lattice sizes; and (iv) discussion of the difference between 
the largest Lyapunov exponents for bipartite and nonbipartite lattices. The computed 
Lyapunov spectra were found to be weakly convex. We have observed no finite 
offset of the smallest positive exponents, and no sign of Goldstonc modes. Both the 
largest Lyapunov exponents and the whole Lyapunov spectra were found to become 
independent of the lattice sizes for sufficiently large lattices. We have given an 
analytical argument explaining this finding. We have found that the largest Lyapunov 
exponent for bipartite and nonbipartite lattices depends differently on the anisotropy 
of the coupling. This is due to the special symmetry of the bipartite lattices with 
respect to the sign change of the two out of three coupling constants. 
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